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ABSTRACT 

Although pulsars are one of the most stable clocks in the universe, many of them are 
observed to 'glitch', i.e. to suddenly increase their spin frequency v with fractional 
increases that range from — s» 10~^^ to 10~^. In this paper we focus on the 'giant' 
glitches, i.e. glitches with fractional increases in the spin rate of the order of — « 10~^, 
that are observed in a sub class of pulsars including the Vela. We show that giant 
glitches can be modelled with a two-fluid hydrodynamical approach. The model is 
based on the formalism for superfluid neutron stars of Andersson and Comer (2006) 
and on the realistic pinning forces of Grill and Pizzochcro (2011). We show that all 
stages of Vela glitches, from the rise to the post-glitch relaxation, can be reproduced 
with a set of physically reasonable parameters and that the sizes and waiting times 
between giant glitches in other pulsars are also consistent with our model. 



1 INTRODUCTION 



The timing of radio pulsars provide us with one of the most 
stable clocks in the universe. Many pulsars exhibit, however, 
sudden increases in the spin rate, known as 'glitches'. To 
date the re have been glitches reported in more than 100 
pulsar^j lEspinoza et al.ll201ll ). with fractional jumps in the 
spin rate of ^ « 10"" to 10"^ 

There is still no clear consensus on the origin of these 
phenomena, but it is thought that glitches (at least the giant 
glitches that are observed in somewhat older, colder pulsars 
such as the Vela) are due to angular momentum being stored 
in a superfluid component of the star that is temporarily 
decoupled from the charged component to which the elec- 
tromagnetic emission is anchored. When the two compo- 
nents recouple there is a sudden transfer of angular momen- 
tum to the crust, whic h gives rise to the observed spinup 
IjAnderson fc Itohlll975l ). A superfluid rotates by forming a 
quantised array of vortices, the distribution of which de- 
termines the rotational profile of the star. The main idea 
at the base of most gl itch models is that vortices can pin 
to the crustal la, t tice llAnderson fc Itohl [19751: lAlpai 1 19771 : 
iPines et al1ll980l : lAJpar et al.lll98ll : [Anderson et all 1198^ . 
allowing for a lag to build up between the superfluid and 
the charged component. 

F ollowing the seminal work of iBavm. Pethick fc PinesI 
|[197l|) several models have been developed to explain the 
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dynamical evolution of the pinned superfluid coupled to the 
crust, mainly to with the intent of explaining the observed 
timescales of days to months in the post-glitch recovery of 
the Vela pulsar. Essentially one can distinguish two classes 
of models, those that assume that the relaxation is due to 
the weak coupling between the superfluid and the crust due 
to the interac tion between free vorti ces and the coulomb lat- 
tice of nuclei l|joneslll990l . 1 1991 [l99Sl and thos e which rely 
on the creep model based on th e seminal work of Alpar et al.l 
(11984) and later adaptations (| Alpar. Langer fc Sauls! Il984 : 
iLink. Epstein fc Bavmlll993l ). In this case the assumption 
is that as a lag develops between the crust and the su- 
perfluid, vortices creep through the crust at a rate that is 
highly temperature dependent, gradually transferring an- 
gular momentum and leading to a steady state spindown. 
However, when the lag approaches a critical value, it is pos- 
sible for an instability to trigger a catastrophic unpinning 
event associated with a sudden transfer of angular momen- 
tum to the crust, after which vortices repin and a new equi- 
librium is reached. These models work remarkably well in 
explaining the post-glitch relaxation, but struggle to ex- 
plain the exact nature of the perturbation which triggers 
it. Such an e vent may be triggered by large temperature 
perturbations iLink fc EpsteinI 1119961 ), caused for exampl e 
by starguakes iBavm fc PinesI (|l97ll ): ICheng et all l| 19921 ). 
the interac tions of the proton vortice s and the crustal mag- 
netic held ISedrakian fc Cordesll 19991 ) or, in the presence of 
strong crustal pinning, a two-stream instability in the su- 
perfluid flow may de yelop, leading to a sudden transf er of 
angular momentum (jGlampedakis fc Andersson! |2009| ) Re- 
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cently, however, IWarszawski fc Melatoa l|201ll ) have shown 
that glitches can be triggered by unpinning avalanches 
of the vortex array, leading to very promising results for 
the overall glitch size and waiting time distributions (see 
also Warszawski fc Melatoa l|2010l '): iMelatos fc Warszawskil 
(l2n09l VL 

An entirely different mechanism was proposed by 
iRudermarJ l|l99ll ), who suggested that vortices pinned to 
the crust stress it to the point of fracture and then move 
outwards with the matter they are pinned t o. Several other 
"starquake" models ha ve been proposed IjEpstein fc LinkI 
I2OOOI : IPranco. Link fc E pstein 2000) but they are all based 
on the main assumption that a starquake triggers the out- 
ward motion of vortices pinned to the matter and the lo- 
cal rearrangement of the crust then causes an increase of 
the angle between the rotation axis and the magnetic axis, 
leading to an increased electromagnetic braking torque on 
the star. Note that starquakes had been invoked as a possi- 
ble explanation for pulsar glitches very soon after the first 
observations, but had been ruled out as they could not ac- 
count for the gia nt glitches of the Vela IjBavm fc Pineal 197ll : 
I RudermarJI 1 9 75^ . This mechanism may, however, play a role 
in other systems and in fact there is some evidence that 
smaller glitches in younger and h otter pulsars such as th e 
Crab may be linked to starquakes (jMiddleditch et al.ll2006l ). 
Furthermore it has been suggested that the glitch may not 
be exclusively linked to vortex dynamics in the crust, but 
that the interaction between the vortices and the supercon- 
ducting flux tubes in the core may contribute to the event 
l|Ruderman. Zhu fc Chenlll998l : ISiderv fc Alpaill2009f) . 

One of the main difficulties in constructing a glitch 
model is, however, the relative lack of realistic calcula- 
tions for the pinning forces. Several calculations exist 
of the pinning energy and of the force per pinning site 
Alpar (1977); Alpar ot al,. (.1984); Epstein_&__Bavm (198g)j_ 
Link fc Epstein (liggih: IPizzochero. Viverit fc Brodia 



lll997l): |joneril99a): iDonati fc Pizzocherol (|2004l l?l l2006h : 
lAvogadro et al.l (|2007l 'l but such calculations of the pinning 
force neglect the finite length of the vortices and rely on 
simplified geometries. Recently, however. Link (2009) has 
analysed the motion of a vortex in a rand o m po tential 
and, even more crucially, iGrill fc Pizzocherd ()201ll ) have 
obtained realistic estimates of the pinning force per unit 
length of the vortex, the quantity that can be directly 
compared to the Magnus force which tends to push the 
vortices out and force the superfluid into corotation with 
the crust. 

IPizzocherd (|2011h has recently shown that a sim- 
ple analytic model, based on the catastrophic unpinning 
paradigm and the realis tic pinnin g forc e s cal culated by 
iGriU fc Pizzocherd l|201lh (see also JGrilll ()201ll )l. can de- 
scribe glitches in the 'Vela' like pulsars, that are typically 
older, exhibit larger glitches (— « 10~^ to lO"'') and al- 
ways s how a decrease in th e frequency derivative after the 
glitch (JEspinoza et al.ll201lh . Essentially the assumption is 
that as the star spins down the vorticity is accumulated in 
the strong pinning region and then released once the result- 
ing Magnus force exceeds the maximum of the pinning force. 

In order to make quantitative predictions for the ob- 
served jump in frequency and subsequent relaxation of the 
crust it is, however, also necessary to describe the evolution 
of the various fluid components that form the NS. In order to 



do this it is convenient to follow a hydrodynamical approach 
that does not deal with vortex motion directly, but rather 
follows the evolution of the separate fluid components. In 
particular we shall make use of the mu ltifiuid formalism de- 
veloped bv lAndersson fc Comerl (120061 ') and recently applied, 
albeit in a mu ch simpler formulation, to the study of pul- 
sar glitches by ISiderv. Passamonti fc AnderssorJ (|2010|). In 
the following we shall thus go beyond the analytic model of 
IPizzocherd (l2011h and show how the realistic pinning results 



of iGrill fc Pizzocherd (|201ll ) can be incorporated into this 
formalism and used to accurately reproduce, for a reason- 
able range of parameters, the observed properties of pulsar 
glitches. Furthermore, as we will discuss in the following sec- 
tions, we will assume that the relaxation is entirely due to 
the mutual friction between the superfluid compone nt and 
the crust. In a sense we are thus taking the view of I Jonea 
( 1992) that the relaxation is due to the dynamics of vortices 
close to corotation with the superfluid, rather than to vortex 
'creep'. Note, however, that the creep model has been shown 
to be consistent with the observed r elaxation timescales of 
the Vela pulsar in the linear regime (| Alpar. Cheng fc Pinea 
1 1989). and we can thus qualitatively account for it by rescal- 
ing the mutual friction parameters, eis will be discussed in 
the following. Finally let us also point out that an effect 
that is neglected in this preliminary model is that of Ekman 
pumping at the crust-core interface, which has been shown, 
together with shear viscosity, to be relevant in explaining the 
post-glitch relaxation times (|van Evsden fc Melatoa l2010f ) 
and should be included in future developments. 

Let us stress that such detailed theoretical work is 
crucial at this time, as the new low frequency radio tele- 
scope LOFAR has just come online an d begun observation s 
of pulsars and fast radio transients (JAlexov et al.l 120111 ). 
As development continues LOFAR, and in the future the 
Square Kilometer Array (SKA), are likely to not only pro- 
vide a wealth of data on both known and new glitch- 
ing pulsars, but also resolve (or at least tightly constrain) 
the glitch rise time, thus allowing us to test our theoret- 
ical understanding of the glitch mechanism ( jSmits et al.l 
2009). Furthermore a glitch may trigger a gravitational wave 
(GW) burst which could be detectable by next generation 
GW observatories llSiderv. Passamonti fc AnderssorJ I2OI0I : 
iBennett. van Evsden fc Melatoaboid ) . By accurately deter- 
mining the glitch time LOFAR and the SKA could be used 
to trigger a directed GW search. 

Finally the constructions of long-term hydrodynamical 
simulations of NS interior could also be extended to model 
not only large pulsar glitches, but more generally pulsar tim- 
ing noise, in order to understand the role that the interior 
dynamics of the star plays compared to magnetospheric pro- 
cesses, which have been shown to be at work in some systems 
IjLvne et ahluOlOf) . This is an issue that is of great impor- 
tance for the current efforts to detects GWs with pulsar 
timing arrays (jHobbs et al.ll20ld V 



2 TWO FLUID EQUATIONS OF MOTION 

We model the pulsar as a two fluid system of superfluid neu- 
trons and a charged component, which consists of the crust 
and the charged particles in the interior. As the neutrons are 
superfluid they will rotate by forming an array of quantised 
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vortices, the interaction of which with the fluids will give 
rise to a weak coupling between the components, known as 
mutual friction. To describe the equations of motion of the 
crust and of the superfluid we shall use the two-fluid formal- 
ism fo r neutron star cores developed bv lAndersson fc Comeij 
(|2006t ). We thus do not consider the equations of motion of 
the vortices, but rather model the macroscopic motion of 
two dynamical degrees of freedom, representing the super- 
fluid neutrons and a neutral conglomerate of protons, elec- 
trons and all non-superfluid components that are strongly 
coupled to them. Assuming that the individual species are 
conserved, we have the standard conservation laws 



ftpx + Vi(Pxt'i) = 



(1) 



where the constituent index x may be either p or n. The 
Euler equations take the form 

where w\^ = nf — uf (y 7^ x), and fi^ — fi^/m^ repre- 
sents the chemical potential (in the following we assume 
that nip — rrin). Moreover, $ represents the gravitational 
potential, and the parameter e^ encodes the entrainment ef- 
fect. The force on the right-hand side of ([2} can be used 
to describe other interacti ons, including dissipative terms 
l|Andersson fc Comerl 120061 ) . We will focus on the vortex- 
mediated mutu al friction force. This m eans that we consider 
a force of form (JAndersson. Siderv fc Comer 2006 1 



/f = KUvPnB'eijki^iw^y + Kn^pnBeijki^ie "^(ifw: 



J klm A,nxy 



(3) 



where Q-' is the angular frequency of the neutron fluid (a hat 
represents a unit vector), k = -^^ represents the quantum 
of circulation and n„ is the vortex number per unit area. 
In particular the quantity nuv is linked to the rotation rate 
of the neutrons and protons by the relation (where f is the 
cylindrical radius): 

d 

fi:n„ = 2 [t7n + en(np - aO] + f^ [S^n + en(!^p - ^n)] (4) 

One can express the mutual friction force in terms 
of a dimensionless "drag" pa rameter TZ such that 
lAndersson. Siderv fc ComejIJOO^ ') 



B 



7^ 



and 



B' ^ 



n" 



(5) 



(6) 



l-f7^2 ' l + 7^2 

and is related to the usual drag parameter 7 by the relation 

Note that for 7?, -C 1 (which is relevant in our context) one 
has S Ri 7^ and B' < B. 

In our model we shall consider the two components to 
be rotating around the same axis defined by Clp. Further- 
more the protons in the crust (to which the magnetic field 
is assumed to be anchored) are bound in ions which form 
a crystalline solid and are connected to the protons in the 
core on an Alfven crossing timescale, which for a neutron 
star core is of the order of a few seconds (or less if the core 
is superconducting) and thus much shorter than the dynam- 
ical timescales we are interested in (except possibly for the 
glitch rise time). Accounting for the elastic forces in the crust 
and Alfven waves in the core is clearly a prohibitive task, so 
as a first approximation we shall account for these effects by 
considering the proton conglomerate to be rigidly rotating. 



We make no such assumptions for the neutrons which will, 
in fact, develop differential rotation if vortices migrate to 
different regions of the star, as can be seen from equation 
©. 

Following ISiderv. Passamonti fc AnderssonI l|2010l ) we 
can write the equations of motion for the angular frequency 
of the two components in the form: 



Q{f) 



Q{r) 



0.n{f) = 



sip(f) 



where we have defined 



Q{f) = pnKn„S(np - Q.n 



-^nl 



and 



A = 



&3 



(7) 
(8) 

(9) 
(10) 



with B the surface magnetic field strength, 6 the inclination 
angle between the field and the rotation axis and /p the 
moment of inertia of the proton fiuid in the star. Note that 
the coupling between the components depends only on the 
dissipative mutual friction coefficient B and not on B' . If we 
now assume rigid rotation for the protons, i.e. a constant 
f2p, we can can multiply equation ([8]) by f^ pn and integrate 
over the volume, thus recasting the system in the form 



tln{f) 



Op 



+ / rQ(:r)-. 



-dV 



Pn 



(11) 



(12) 



3 THE PINNING FORCE 

Although vortex pinning in the crust has been suggested 
as the main mecha nism behind pulsar glitc hes more than 
twenty years ago bv I Anderson fc Itohl fl975), until recently 
very little work has been devoted to obtaining realistic esti- 
mates of the pinning force. The main difficulty lies in the 
fact that, although there have been several estimates of 
the pinning energ y and thus of the force per pinning site 
(Alp ar et al.lll984 ). it is in fact the pinning force per unit 
length acting on a vortex that is the relevant quantity to 
compare with the Magn us force if one is to und erstand when 
a vortex line can unpin IjAnderson et al.l 19821). 

In fact it has been argued by Ijoned ([1997|) that if 
one considers an infinitely long vortex line and averages 
over the various orientations with respect to the lattice, 
the pinning force will be negligible. Recent calculations by 
iGrill fc Pizzocherd l|201H ) have shown that for realistic con- 
figurations in which one considers a vortex to be rigid over 
lengthscales of 100-1000 Wigner-Seitz radii, the averaging 
process over mesoscopic lengthscales and different orienta- 
tions of the lattice does indeed reduce the strength of the 
pinning force with respect to previous estimates based on 
the pinning force per pinning site. However, the values ob- 
tained in these calculations, that are in fact the first realistic 
estimate of the pinning force per unit length, are s till large 
enough to explain pulsar glitches (|Pizzochercil l201lh . 

It is beyond the scope of the present work to 
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make a comparison between various superfluid gap mod- 
els and equations of state (howeve r see the discussion in 
IPizzochero. Seveso fc Haskelll l|201ll '). who find that the main 
quahtative features of the model remain unchanged), as at 
present we are mainly interested in understanding how to re- 
produce the main features of a glitch with a simplified large 
scale hydrodynamical neutron star model. Furthermore we 
shall use a Newtonian model to describe the mutual friction 
dissipation between the two components, as at present the 
relativistic prescription for doing this is not fully developed. 
Given the inaccuracies inherent to a Newtonian calculation 
of a neutron star model any comparison between realistic 
equations of state would be at least dubious. Furthermore 
there are, as we shall see, huge uncertainties associated with 
the values of the mutual friction coupling parameters that 
would make such a calculation meaningless. 

In this work we shall focus on the case of continu- 
ous v ortices that thread the whole star, as in (|Pizzocherd 
2011^■ and extract the main features of the calculation of 
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Pizzochero. Seveso &: Haskelll l|201ll ') to construct a simpli- 



fied model for the maximum lag that can be developed. In 
particular, as the maximum pinning force is obtained by bal- 
ancing the Magnus force over a whole vortex the maximum 
lag between the two components (AQ) will have cylindri- 
cal symmetry and be a function of the cylindrical radius 
f = rsm{6) only. Close to the rotation axis the radial be- 
haviour is thus due to the Magnus force, as the vortices 
encounter a roughly similar number of pinning sites. The 
maximum lag thus falls off as « 1/f until it reaches a mini- 
mum value of An « 10~* rad/ s at the base of the crust (note 
that i n the analytic model of IPizzochero. Seveso fc Haskelll 
l|201ll ) the minimum is at somewhat higher densities, deeper 
in the outer core). The maximum lag then rises steeply in 
the equatorial region as the vortices are now completely con- 
tained in the pinning region and reaches a maximum value 
of An ~ 10~^ rad/s, in the inner crust. We model this be- 
haviour as a linear rise of the lag between the minimum 
and maximum values. As we continue to move outward the 
lag then decreases linearly again as the vortices annihilate 
in the outer regions of the crust where the neutrons are no 
longer superfluid. A schematic representation of this is given 
in figure [l] 
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Figure 1. The top panel shows our approximation for the critical 
unpinning lag for a 12 Km, 1.4 Mq star, modelled as an n = 1 
Newtonian polytrope. In the bottom panel we show for compar- 
ison the critical unpinning lag for a realistic model of a 12 km , 
1.4 M0 using the SLy equation of state (C haban at et al.lll998l) 
as discussed in IPizzochero. Seveso fc Haskelll l|201lf l. 



4 EQUATION OF STATE AND DRAG 
PARAMETERS 

Let us discuss the fiducial NS model that we adopt for our 
calculations. We consider a 1.4 Mq, with a radius _R = 12 
km, and assume that the base of the pinning region is at r = 
11.15 km, at a density p = 9. 6 x 10^^, which is rougly con- 
sisten t with the estimates of IPizzochero. Seveso fc Haskelll 
l|201ll ). We take the background equation of state of our 
neutron star to be a n = 1 polytrope, which allows us to 
make use of the explicit solution for the density profile 



P = Pc 
with 



Txr , 

a = -— - and 
R 



Pc = 



ttM 
4R3 



(13) 



(14) 



This gives us a value It = 1.57 x 10*^ for the total moment of 
inertia. However we need to consid er the proton and neutron 
densit ies separately, so, following iReisenegger fc GoldreichI 
(Il992f ). we shall consider a proton fraction Xp of the form: 



= ^= 0.05 



10i4g/cm^ 



(15) 



where p = pp + pn is the total density. 

Finally in a fully consistent model one should be 
able to obtain the entrainment coefficients directly from 
the equation of state. There are however few realistic 
fully consistent equations of state for superfluid neu- 
tron star matter. Entrainment coefficients have, however, 
been obtained for the cores of sta rs containing hyperons 
(|Gusakov. Kantor fc Haenselll2009al la) and for the c rust one 
can r e ly on calculations of the proton effective mass JChameJ 
l2006l : IChamel fc Carterl l2006l : IChamel. Pearson fc Gorielvl 



Modelling pulsar glitches 5 



le-03 FT 



le-04 - 



■a le-06 ? 



^ le-07 - 

& 

Q 



le-09 r 



le-10 - 



I I I I I I I I I I I I I I I I I I I 
'■" "O 1 2 3 4 5 6 7 8 9 10 11 12 
Cylindrical Radius (km) 

Figure 2. The average drag parameter 7^ as a function of cylin- 
drical radius, obtained by averaging the core contribution TZc and 
crust contribution TZp over the length of a vortex. Note that we 
consider only the drag due to phonon excitation (TZp) in the crust 
as we are assuming that Kelvin waves can only be excited once 
the critical unpinning lag has been exceeded. The horizontal line 
corresponds to TZ of the order of 10~®, which gives a coupling 
timescale of the order of 1 minute. 



I2OI0D which is related to our e ntrainment coefficients by 
tlie relation (|Prix. Comer fc Andersson.2002 ): 



(16) 



where nip is the bare proton mass and m* is the proton 
effective mass. Recent calculations Chamel ^2006) suggest 
that the proton effective mass is slightly lower than the bare 
mass in the core but can be larger in the crust. This means 
that the entrain ment parame t ers w ill vanish close to the 
base of the crust ICarter et al.l (12006! ) . As this is the region 
we are mainly interested in for our evolutions and given the 
qualitative nature of this work, we shall take Ep = e„ = 0. 

The discussion of the mutual friction coefficients is more 
complex. The mechanisms that give rise to mutual fric- 
tion are in fact different in different regions of the star. 
In the core we expect electron scattering of vortices to 
be at work, coupling neutron s and protons on a timescale 
of T gj 10 rotational periods (lAlpar. Langer fc Sauls! 1 1984 
lAndersson. Siderv fc Conieill2006f ) (which for the Vela pul- 
sar would give a timescale of slightly less than a second). 
Another possibility is that, if the protons are in a type 
II superconducting state, the interaction between fluxtubes 
and vortices will c ouple the two compo nents on an even 
shorter timescale f|Siderv fc Alparl 120091 ) . although if the 
flux tube tangle is sufficiently strong to "pin" the vor- 
tices this may lead to the opposite being true and the 
core being effectively decoupled for most o f the evolution 
l|Ruderman. Zhu fc Chenlll998l : I Linkll2003l ). We shaU dis- 
cuss this possibility in more detail in the following. 

In the crust, on t he other hand, it has been shown b y 

iJonej (|1992D (see also (|Joneslll990l : lEpstein fc Bavmlll992l )) 
that coupling to sound waves in the lattice will be the main 
source of dissipation for low velocities of the vortices with 



Table 1 ■ Drag parameters in the crustal regions defined by 
IINegele fc V authcrin 1973,). In the last two lines we have applied 
a reduction factor 5 = 10""*, as described in the text. 7?.p indi- 
cates the drag parameter due to phonon excitations and TZ^ that 
due to Kelvon excitations. 

Zone 12 3 4 5 



p/lO^^'g/cmS 


0.015 


0.096 


0.34 


0.78 


1.3 


TZp/W-^ 


1.6 


0.7 


5.8 


0.025 


0.0 


Tlk/lO^ 


290 


5.8 


340 


0.085 


0.0 


Sp/10-5 


1.6 


0.7 


5.8 


0.025 


0.0 


8fc/10-2 


0.003 


0.18 


0.003 


11.64 


0.0 


5 = 10-4 



Bp/10-^ 



1.6 
15.5 



0.7 
2.9 



5.8 
13.5 



0.025 
0.04 



0.0 
0.0 



respect to the lattice (^ 10^ cm/s) while dissipation due to 
the excitation of Kelvin waves in the vortices will dominate 
for large relative velocities. 

In table [1] we show the values of the drag parame ters in 
the cr ust, for both these interactions, obtained from Ijoned 
(|l990l ) by using recent values of t he vortex-nucleus inter- 
action IjDonati fc Pizzocherdl2006l ). Unfortunately there is 
a large uncertainty associated with this calculation, as one 
must not only evaluate the energy dissipated by the inter- 
action of the vortex with a single pinning site, but also sum 
coherentl y over the lengthscale asso ciate d with vorte x rigid- 
ity fwhich lGriU fc Pizzocherol l|201ll ) and lCrilll l|201ll ) find to 
be of the order of 10^ — 10^ Wigner-Seitz radii). The summa- 
tion proc edure leads to a poorly constrained reduction factor 
S, which (|joneslll990l ) finds to be of the order S » 10"^ Thi s 
is consistent with the results of iGrill fc Pizzocherol l|201l[ ). 
who find that averaging over mesoscopic vortex length-scales 
and different crystal orientations leads to a reduction of up 
to two orders of magnitude in the pinning force per unit 
length, with respect from previous estimates that were ob- 
tained directly from the pinning force per pinning site. Given 
the scalings in the expressions for the drag parameters of 
Ijonesi (|l990 n. this would lead to a reduction factor 5 ~ 10"* 
for the drag parameters. There is thus a large uncertainty 
associated with the values in tableland, furthermore, drag 
parameters 7?. >> 1 are dubious as f or such a str ong in- 
teracti on the perturbative trea tment of lJoned (Il99d ) breaks 
down (|Epstein fc Bavmlll992l '). This suggests that for our 
qualitative study we may take a constant drag parameter 
throughout the crust and study how the results depend on 
its variation. 

In particular we shall consider values in the region 
Bp ~ 10^^" for the weak drag due to the interaction 
with lattice phonons and Bk ~ 10"'^ for the strong drag 
due to Kelvin wave excitation. For the core, on the other 
hand we shall consider electron sca ttering off vortex cores 
as th e main source of dissipation IjAlpar. Langer fc Saulj 
1 19841 ). In this case we take for the mu tual friction coeffi- 
cient (|Andersson. Siderv fc Comerll2006r ): 



Bc = Ax 10" 






(17) 



and 



where m* is the effective proton mass, Jm* = m* - 

pi4 is the density in units of 10^* g/cm^. However, given that 
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we are already taking the drag parameter to be a constant 
in the crust we shall also consider constant drag in the core. 
This choice is more than adequate given that we are not 
using a realistic prescription for the density profile or the ef- 
fective mass. We shall thus consider drag parameters in the 



region of Be ~ TZc 



5 X 10" 



Note that microscopically 



the exact nature of the superc onductor in the neutron star 
interior is still open to debate ( Jonos|[2006|). The values we 
use are consist ent with the cor e being in a typ e I supercon - 
ducting state (ISedrakianll2005l ') (although see l|jonej 120061 ') 
for a discussion of strong drag in a type I superconductor) 
and with the possibility that one has a type II superconduc- 
tor and the drag is stro ng, either as a c onsequence of the 
magnetic field geometry (jSiderv fc Alpanl2009l ') or due to a 
significant number of vortices cutting though the magnetic 
flux tubes. We do not consider here the possibility that vor- 
tices are strongly pinned to flux tubes in a type II supercon- 
ductor and the most of the core is thus decoupled from the 
crust. This is clearly an interesting possibility which shall 
be the object of future work. 

In our model the exact value of the drag parameter 
acting on a vortex section depends critically on the region 
we are considering, as not only do the parameters depend 
on density, but the very nature of the drag varies from crust 
to core. We shall therefore integrate over a whole vortex 
and deflne an "effective" drag coefficient which depends only 
on cylindrical radius (as shown in figure [2]) by averaging 
over the vortex length. Our effective drag parameter thus 
depends mainly on how sizable a portion of the vortex is 
immersed in the core (i.e. the stronger drag region) rather 
than in the crust. 



5 DESCRIPTION OF THE MODEL 

5.1 Core 

Let us now describe the evolution of the system as it spins 
down due to the electromagnetic torque. First of all let us 
focus on the central regions of the star. Here the vortices 
will stretch through the fluid core and only a small portion 
will experience pinning to the crustal lattice. This setup 
is in fact very similar to that used to study the spin-up 
of superfluid helium in a container. In this case the vor- 
tices are only pinned to the surface of the container and 
unpin once the Magnus force integrated over a vortex is 
approximately equal in magn i tude t o the tension (see e.g. 
lAdams. Cieolak fc GlabersonI (|l984h : ISoninI ()l987l )): 



Tl 



■inf- 



47r 



(18) 



with b the intervortex spacing and f the vortex core ra- 
dius. This is in fact natural as we can expect that, once 
the Magnus force exceeds the tension the vortex will no 
longer behave rigidly but we can expect reconnection and 
Kelvin waves to be excited, possibly developing turbu- 
lence and leading to a number of vo rtices to unpin and 
begin to "creep" radially outwards (jSchwarj [l98J) (al- 
though bundles of vortices could be significantly more rigid 
iJRuderman fc Sutherland! 1 19731 )). For average parameters 
the Magnus force in the core will exceed the line tension 
almost immediately once a lag begins to develop (Af2c ~ 



10~^^), leading to the conclusion that (as in ''He) unpinned 
vorticity dominates for most of the time, and re-pinning can 
only take place if the the two components are essentially co- 
moving (i.e. AQ, < AQc)- Note the the situation is radically 
different in the crust, where the vortex line is fully immersed 
in the lattice and is pinned in several points. In this case we 
can expect it to move freely only once the maximum pinning 
force of section [3] is exceeded. 

Let us, however, stress that the assumption of free vor- 
tices is not crucial for our model. One can, in fact, easily ac- 
count for the fact that possibly not all vortices will be free 
to move, but at any given time only a small f raction will 
be free, and subject to the drag force. Following 'Jahan-Miril 
(|2005l . l200m we shall assume that the instantaneous number 
density of free vortices rim is given 



Tim — ^Tly 



(19) 



where n„ is the total number density of vortices, and ^ is the 
fraction of unpinned vortices at a given time. By now aver- 
aging over time we can obtain an effective drag parameter 
Be = ^B. 

For the standard "creep " model (lAlpar et al.l 11984 : 
I Link. Epstein fc Bavml 119921 ) the unpinning probability 
takes the form 



C«exp --^ 1 



(20) 



where Ep is the Energy barrier that the vortex needs to over- 
come, and T is temperature, lu is the lag between the com- 
ponents, and ujc is the critical lag. Unfortunately there are 
huge theoretical uncertainties on the unpinning probability 
and observational constraints are also very model depen- 
dent. However by varying ^ and thus reducing the 'effective' 
drag parameter we can account for vortex "creep" in the lin- 
ear regime, which has been used to interpret the re laxation 
phase of Vela glitches (jAlpar. Cheng fc Pineg|l989l ). 

The details of the steady state are, however, not cru- 
cial for our formulation, which is fortunate as not only vor- 
tex creep but also vortex repinning is still poorly under- 
stoo d, although efforts are being m ade to tackle this prob- 
lem (|Sedrakianlll995l : I Linl3l2009l '). Eventually the system 
will settle down to an equilibrium state in which both fiu- 
ids are spinning down at the same rate and the lag can be 
estimated, from the equations in ^ and ((8|, to be 



An{f) 



n 



2B(r)n 



(21) 



where B is the effective (averaged) mutual friction param- 
eter and we have neglected the effect of entrainment and 
differential rotation. We will see in the following that this is 
actually a very good approximation to the equilibrium lag. 



5.2 Strong pinning region 

The repinning region, i.e. the equatorial region in which the 
pinning force rises steeply, deserves a separate discussion. 
As the vortices encounter what is, effectively, a barrier, we 
shall assume that they repin almost immediately. As we dis- 
cuss below this is effectively a situation in which the vortices 
'creep' steadily outwards, i.e. move outwards with, on aver- 
age, a low velocity. This suggests that Kelvin-phonon inter- 
actions will be strongly suppressed and that we can con- 



Modelling pulsar glitches 7 



sider weak mutual friction coefficients, due to interactions 
with lattice phonons, in the region of Bp ^ TZp ^ 10"^'', as 
described in the previous sections. 

In f act, as the d rag description is not valid close to re- 
pinning l) Linkll2009l l. one could argue that the form of the 
mutual friction in (JS)) is also not valid. Strictly speaking 
this is true but, as the full problem of the motion of a vor- 
tex in a strong pinning potential is at present intractable, 
we shall assume that the coupling given by mutual friction 
with the weak drag given by the interaction with phonons 
in the lattice is a good approximation. The whole discussion 
could also be made rigorous by followi ng the approach of 
lAndersson. Haskell fc SamuelssorJ l|2011f ) and assuming that 
an extra 'pinning' force is acting on the vortices. This leads 
to an effective mutual friction parameter 



1 + 7Z2 



7^ + 



ail + a±TZ 



(22) 



which now depends on the relative velocity w and the com- 
ponents of the effective force that describes the motion 
through the pinning potential, ay (along the relative flow) 
and a± (perpendicular to the relative flow) . Note that this is 
not the actual "pinning" force, but rather an effective force 
that mimics the effect of creep. The vortex velocities then 
take the form 



W|| = ayT?, — w — a± 

and 

v± = ay + {w + a±) Tl 

and in particular one finds that 



B: 






(23) 



(24) 



(25) 



However, given that we do not have a consistent descrip- 
tion of an effective pinning force (but only the maximum 
pinning force) we shall simply assume that all vortices are 
pinned below the critical velocity, and then relax to their 
steady state configuration once the critical unpinning lag 
is exceeded. This is equivalent to assuming that below the 
critical velocity only a very small fraction of the vortices 
can 'creep', which is reasonable due to the sharp rise in the 
pinning force. 

This leads to the neutrons relaxing to their steady state 
rotational rate on a timescale Tr « l/2ftB ~ 10* s (for B ~ 
10~^°), during which the protons have spun down and the 
unpinning "front" has moved out along the strong pinning 
equatorial region. Essentially the unpinning front moves ra- 
dially out at a speed Vr — Ar/rg with Tg = Aflmax/^, which 
leads to a speed of approximately Vr ~ tlAr/AQmax ~ 10~^ 
cm/s, where Ar = 400m is the width of the strong pinning 
region, and AQrnax ~ 10~^ is the maximum in the lag in the 
equatorial region. This means that the vortices that have un- 
pinned from parts of the star closer to the rotation axis are 
clustered in this region, over a lengthscale I « TrVr ~ 10* 
cm. We shall assume that it is this region that gives rise to 
the glitch once the maximum lag is exceeded and the vor- 
tices can no longer be pinned, thus no longer 'creep', but are 
free to escape radially outwards. 



5.3 Glitch 

Once the maximum lag has been reached we shall assume 
that the pinning force can no longer contrast the Magnus 
force and that the vortices move out. They very rapidly 
reach a state of co-rotation with the neutrons, leading to 
velocities of « 10* cm/s with respect to the lattice. In 
this regime Kelvin-phonon dissipation will certainly dom- 
inate and recouple the neutrons and protons on a very 
short timescale, giving rise to the rapid exchange of angular 
momentum that characterises the glitch. Consistently with 
our prescription for the electron-phonon dissipation we shall 
take the mutual friction coefficient for Kelvin-phonons to be 
Bk ~ 10~^, and assume that acts only in the region that has 
not relaxed to a steady state, i.e. the region in which the lag 
is greater than the equilibrium value: 



An > 



n 



2Bpn 



(26) 



We can estimate the relaxation timescale in the crust as t = 
l/2QBp and assume that the pinning front moves through 
the strong pinning region over a period of 3 years, i.e. with 
a velocity « ~ 4 x 10~ cm/s. The region that has not yet 
relaxed will then have a thickness of 



3 X 10"® 



(27) 



For a drag parameter Bp ~ 10^^" one would then have ap- 
proximately all of the strong pinning region involved in the 
glitch. A more quantitative analysis is shown in figure [3] 
where we show the edge of the non relaxed region as a func- 
tion of the drag parameter, obtained by measuring where 
the lag between neutrons and protons deviates from the 
analytic estimate in (|26|) . We cannot extend the numerical 
analysis to very weak drag, nevertheless we see that the re- 
gion becomes larger for weaker drag and that extrapolating 
our results we would have the whole crust involved in the 
glitch for Bp ~ 2 X 10"^". The extrapolation was obtained 
by fitting a function of the form, inspired by the estimate 
in (|26|l . f{x) — a — b/Bp, where the best fit parameters 
were a = 2330.45 ± 26.14 m and 6 = 1911.67 ± 203.9 m, 
which are in reasonable agreement (within factors of a few) 
with the analytical estimate. We shall thus assume that for 
Bp < 2 X 10"^" Kelvin oscillations of the vortices can be 
excited in the whole strong pinning region. 



5.4 Post-glitch relaxation 

Subsequent to the glitch neutrons and protons will still be 
in approximate co-rotation only in the interior regions that 
are coupled on timescales shorter than the glitch rise time. 
In the exterior core the lag will have decreased by a factor 
AQ ~ 10"* due to the glitch itself. The neutrons and pro- 
tons then return to their steady state configurations with 
the timescale appropriate to the mutual friction parameter 
B at the cylindrical radius we are considering, thus naturally 
recoupling the core on timescales that range from the glitch 
rise time in the innermost regions to months in the outer- 
most regions (i.e. f ft; 11 km). As far as the strong pinning 
region is concerned we shall assume that the vortices repin 
after the glitch, and gradually unpin (or begin to 'creep') 
again as a lag builds up due to the protons spinning down. 
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Figure 3. We plot the extent of the region in which Kelvin os- 
cillations can be excited as a function of the weak drag param- 
eter TZ-p. Note that the whole strong pinning region is 400 me- 
ters thick, as indicated by the horizontal line. We extrapolate 
the value of the drag parameter for which the whole region is in- 
volved by fitting a function of the form, inspired by the estimate 
in II26I I. f{x) = a + b/{'R,p/10~^^), where the best fit parameters 
were a = 34.77 it 13.7 m and b = 955.8 it 101.9 m. As we can 
see the whole strong pinning region participates in the glitch for 
TZp < 2 X 10""^'', as expected from the analytical estimates. 



As already discussed the details of the repinning are highly 
uncertain, but given the weak coupling in the crust the de- 
tails of such a mechanism will only impact on the longer 
relaxation timescales. 



6 RESULTS 

In order to set up the pre-glitch conditions for our system 
we impose that the two fluids rotate at the same rate and 
then begin to evolve in time. This would clearly not be the 
correct condition if, say, we are considering the situation 
immediately following a previous glitch. However, as we are 
setting up the situation for the next glitch we are only in- 
terested in obtaining a steady state on long timescales (i.e. 
the inter-glitch timescale, which is approximately 3 years 
for the Vela pulsar). As we intend to evolve the system for 
approximately 3 years we cannot evolve the whole interior, 
which would require us to deal with coupling timescales of 
seconds, but assume that it is always coupled to the crust 
by including its moment of inertia in that of the charged 
component, and only evolve the region between f = 11.14 
Km and r = 11.55 Km (i.e. the outer core in which the value 
of the drag parameter becomes small and varies rapidly, and 
the strong pinning region of the inner crust). We neglect the 
outer crust, which we do not expect to be a bad approxi- 
mation as its moment of inertia is small compared to that 
of the other regions of the star. The equations of motion 
are thus integrated only up to the maximum of the pinning 
force, where we impose that the spatial derivative of the 
neutron angular velocity must vanish, i.e. ^^ = 0, which 
is appropriate if the vortices are still pinned at that point. 
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Figure 4. We compare, in the strong pinning region of the crust, 
the analytical prediction for the lag Af2 between protons and 
neutrons to the value obtained in our simulations just before the 
glitch. The calculation was performed for Hp = 10~^. As can be 
seen from the figure, the analytical estimate is a very good ap- 
proximation, except for the outer region in which the system has 
not yet reached equilibrium. It will thus be a good approximation 
for the core, in which the coupling timescale between neutrons and 
protons is considerably shorter than in the crust. 



At every time-step we thus evaluate the integral in (lll|l and 
then evolve the coupled equations in (|lllll2p forward in time 
with a four step Runge-Kutta algorithm. 

Once the maximum of the lag in the strong pinning re- 
gion is reached we initiate the glitch by switching to strong 
Kelvon drag (Bk) in the region that has not relaxed, i.e. the 
region for which Af2 > 2B n ' ^^ ^'^^^ point we track the 
evolution of the whole core by extending the numerical grid 
and filling the region that had not evolved with the prescrip- 
tion Afi = 



2Bcn 



As can be seen in figure (|4} this is indeed 
a good approximation for the initial condition, given that 
the timescale on which the interior reaches a steady state is 
considerably shorter than the inter-glitch timescale. When 
the strong drag regions have reached their equilibrium lag 
we switch off the Kelvon drag and assume that the vortices 
have repinned. We then follow the immediate post-glitch re- 
laxation by evolving the whole star over a timescale of days. 
To summarise, we set up a system of co-rotating flu- 
ids with pinned vortices in the crust and evolve it in time. 
As the lag increases the depinning front moves thorugh the 
crust. Once the maximum critical lag is reached the vortices 
move rapidly outwards and strong Kelvon mutual friction re- 
couples the components, giving rise to the glitch. After the 
glitch the lag has increased throughout the star and each 
region relaxes to equilibrium on a timescale determined by 
the local values of the (average) mutual friction. 



6.1 The Vela Pulsar 

Let us flrst of all examine the case of the Vela Pulsar, 
which is the prototype system for giant glitches. The Vela 
(PSR B0833-45 or PSR J0835-4510) has a spin frequency 
V ~ 11.19 Hz and spin-down rate v ~ —1.55 x 10~^^ 
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Figure 5. In the top panel we plot the size of the glitch (at the 
end of the rise) and in the bottom panel the rise time, both as 
a function of the strong Kelvin drag parameter 7?.^. For these 



simulations we have take 7?.„ 



10- 



and T^c = 5 X 10" 



Hz s~ . Giant glitches are observed roughly every two to 
three years, and have relative frequency jumps of the or- 
der ^vjv ~ 10^^. The spin- up is instantaneous to the ac- 
curacy of the data, with the best constraint being an up- 
per li mit of 40s for the ri se time obtained from the 2000 
glitch (|Dodson et al.1 I200H ) (a similar upper limit of 30 s 
was obtained from the analysis of the 2004 g litch, but was 
less s ignificant due to the quality of the data (jDodson et alj 
[20071)). The glitch is usually fitted to a model consisting 
of permanent steps in the frequency and frequency deriva- 
tive and a series of transient terms that decay exponen- 
tially. It is well known that at least three are required to 
fit the da ta, with decay t imescales that range from months 
to hours (|Flanaganl Il996l ) . Recent observations of the 2000 
and 2004 glitch have shown that an additional term is re- 
quired on short timescales, with a decay time of approx- 
imately a minute. The fitted values for the relaxation of 
these two glitches are shown in table (2] The spindown rate 
always increases after a glitch with larger relative increases 



Table 2. The fitted valu es for the relaxation o f the Vela 2000 and 
Vela 2004 glitches, from iDodson et"all l l200ll ) and iDodson et"all 
|2007). After removing the pre-glitch spindown the fit on the 
residuals is performed with a function of the type /(t) = Ap_F + 
ApFt-H^^/,exp(-t/rO. 



2000 



2004 



ApF/Hz 


3.45435E-05 


2.2865E-05 


ApF/Hz s-i 


-1.0482E-13 


-1.0326E-13 


/i/lO-^Hz 


0.02 


54 


/2/IO-6HZ 


0.31 


0.21 


/3/IO-6HZ 


0.193 


0.13 


/4/10-SHz 


0.2362 


0.16 


n 


1.2±0.2 mins 


1±0.2 mins 


T1 


0.53 days 


0.23 days 


T3 


3.29 days 


2.10 days 


T4 


19.07 days 


26.14 days 



(up to a factor of 10) which decay on the shorter timescales 
and smaller (a few percent) increases that decay on the 
longer timescales of days and months. Note that a further 
two glitches were detected in 2008 and 2010, but no timing 
analysis has yet been pubHshed. 

Co nsistently with the results of Grill fc Pizzocherd 

(|201ll ). lGrilil (|201ll ') and (|Pizzochercll201ll ) we take the max- 
imum of the critical unpinning lag to be Ailc = 10"'^, which 
naturally leads to a glitch recurrence time of roughly 3 years 
(the exact time depends on the choice of drag parameters, 
but is always approximately 1100 days). Our simulations can 
then reproduce glitches with fractional rises ISv jv ~ 10~^ 
for a range of parameters. In figure [5] we show the rise time 
and the 'instantaneous' size of the glitch, i.e the size at the 
end of the rise, as a function of the strong Kelvin drag pa- 
rameter 7?.strong. As expected the larger the drag parameter 
the shorter the rise time and the larger the glitch. This is 
simply due to the fact that the faster the glitch (or rather 
the stronger the Kelvin drag parameter with respect to the 
drag in the core) the smaller the fraction of the core su- 
perfluid neutrons that can remain coupled to the crust and 
contribute to its moment of inertia during the glitch. 

However a strong Kelvin drag will lead to rise times 
of the order of seconds, well below the observational upper 
limit of 40 s and thus well below the observational capabili- 
ties of current radio telescopes that can, for the Vela pulsar, 
at best deal with 10 s folds. The rise will then be followed 
by a rapid relaxation. In fact if we then extract the glitch 
size 40 s after the glitch is initiated, as shown in figure[6l we 
see that the situation is reversed. The stronger the drag the 
smaller the glitch after 40 s, as an initially larger jump in fre- 
quency relaxes faster and leads, in fact, to a lower frequency 
after 40 s, which is illustrated in figure[7l In figure[6]we also 
see that stronger phonon drag parameters in the core gives 
rise to smaller glitches, as more of the crust has relaxed to 
its equilibrium configuration and less angular momentum is 
available to be exchanged. 

While the magnitude of the strong Kelvin drag affects 
mainly the glitch size and rise time, the strength of the drag 
in the core will affect, as shown in figure [8] and [9] the size of 
the glitch itself (as it determines the amount of core super- 
fluid that participates in the glitch) but also crucially affect 
the short term post glitch relaxation. In figure [9] we can see 
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Figure 6. The size of tlie glitcli, extracted 40 s after it is initiated, 
as a function of tlie drag parameters. In the top panel we plot the 
glitch size as a function of the strong drag parameter T^j; . In this 
case the stronger the drag the smaller the glitch, as an initially 
larger glitch relaxes faster and is, in fact, smaller after 40 s. The 
remaining drag parameters were taken to be TZp = 5 X 10~^^ and 
7^c = 4 X 10"*. In the bottom panel we show the glitch size as a 
function of the weak drag parameter 'R.p . Again for larger values 
of the drag we obtain a smaller glitch, as in this case more of 
the fluid in the crust has relaxed to its equilibrium configuration 
prior to the glitch and less angular momentum is available to 
be exchanged. The remaining drag parameters were taken to be 
Tefe = 6 X 10"" and ■Rc=iy. 10"'*. 
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Figure 7. The first 40 s of the glitch for two values of the 
Kelvin drag parameter TZu- A stronger drag (i.e. a shorter cou- 
pling timescale) gives rise to an initially larger glitch that however 
rapidly relaxes to a lower spin rate than that of a glitch involving 
a weaker drag parameter. 
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Figure 8. The size of the glitch, extracted 40 s after it is initiated, 
as a function of the core drag parameter. Weaker drag parameters 
give rise to larger glitches as less of the core can contribute to the 
moment of inertia during a glitch. The remaining drag parameters 
were set at T^p = 5 X 10"^^ and Uk = 1 x 10~^. 



that a weaker drag in the core leads to less of the superfluid 
coupling to the crust and thus to a larger glitch, that how- 
ever relaxes, on a timescale of several minutes, more than 
in the case of stronger coupling. Unfortunately the short 
timescale component of the relaxation is not strongly con- 
strained by observations, as it was only measurable in the 
2000 and 2004 glitches, with vastly different magnitudes, 
as can be seen from figures [TT] and 1101 The results in fig- 
ure [TT] would however indicate that a core drag parameter 
7^ ~ 10^^ is still marginally consistent with observation, 
but that there is no evidence for a much weaker drag, due 



for example to the fact that most vortices are pinned (either 
to the crust or to fiux tubes in the core) and only a small 
fraction of them can creep. It would thus appear that the 
observations of the short term relaxation are consistent with 
a mutual friction drag in the range R « 10~* — lO"'^, con- 
sistent with theoretical expectations for electron scattering 
off vortex cores. This is encouraging, given the approximate 
treatment of the drag parameter in this work. Note, how- 
ever, that quantitative statements are not easy to make as 
there is a large degeneracy between the various parameters 
that enter the model. For example in figure[T2]we show that 
the uncertainty on the time of the glitch, which is of about 
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Figure 9. We plot the post-glitch residuals after subtracting the 
pre-glitch spindown. The top panel shows the first f5 minutes of 
the relaxation for three values of the core drag parameter T^corc • 
We can see that a weaker drag in the core leads to less of the 
superfluid coupling to the crust and thus to a larger glitch, that 
however relaxes, on a timescale of several minutes, more than in 
the case of stronger coupling. 
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Figure 10. The post-glitch residuals, rescaled as in figurc lTTI but 
for a longer timescale of 2 days. 



30 s, can lead to an uncertainty on the value of the drag 
parameters. 

In figure [10] we show the longer timescale components 
of the relaxation. These also appear consistent with a strong 
core drag due to electron scattering, although a slower spin- 
down may be preferred. Quantitative statements regarding 
the longer relaxation timescales (days to months) are, how- 
ever, hindered by computational and theoretical issues. On 
the one hand evolving the full system of coupled equations 
for longer than a few days is computationally challenging, 
on the other theoretical uncertainties regarding repinning 
after a glitch and creep, begin to have a significant impact 
on the results for these longer timescales, while they do not 
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Figure 11. We plot the post-glitch residuals after subtracting 
the pre-glitch spin-down and having scaled all the results so that 
the glitch has decayed to the same frequency after 10 minutes. 
The bottom panel is a zoom in of the top panel. It is clear that 
stronger drag parameters provide a better fit to the Vela 2000 
relaxation fit (for which the data was of better quality) and that 
very weak drag parameters are still excluded by the Vela 2004 
relaxation fit. It would appear that the observations are consistent 
with a drag parameter in the range R fs 10"'* — 10"'^, as expected 
theoretically for electron scattering off vortex cores. A very weak 
drag, due to only a small fraction of the vortices 'creeping', would 
not appear to be consistent with the short timescale components 
of the relaxation. 



affect short timescales of minutes or hours. In order to ob- 
tain truly quantitative results for longer timescales it would 
be necessary to address these theoretical issues and also to 
include a realistic density dependence of the drag param- 
eters. Finally other effects are likely to play a role in the 
relaxati on timescale, notably frictio n at the crust-core in- 
terface (|van Evsden fc Melatosll20ld ). the inclusion of which 
would require us to relax the rigid rotation assumption for 
the charged component, which is beyond the scope of this 
paper but will be the focus of future work. 
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Figure 12. We plot the post-glitch residuals, rescaled as in figure 
nil to illustrate the degeneracy between the assumed glitch epoch 
and the strength of the drag in the core. We show the effect of 
assuming that the glitch has occurred 30 s before the assumed 
epoch. We can see that the curve for TZc = 4 X lO"* now produces 
a smaller glitch and gives a better fit to the Vela 2000 data. The 
curve thus appears closer to that given by taking TZc = 9 X 10~^ 
but not adjusting the glitch epoch. 



7 GIANT GLITCHERS 

So far we have only applied our model to the giant glitches of 
the Vela pulsar. Giant glitches have, however, been observed 
in several other pulsars and we would expect them to also 
occur when the system reaches the maximum critical lag 
for unpinning. As already mentioned this needs not be the 
case for smaller glitches that are likely to be due to random 
unpinning and may be, for example, described in terms o f 
vortex avalanche dynamics ijWarszawski fc Melatod l2011f l . 
In fact recent observations support the view that there is 
a bimodal distribution in the glitch size, with two different 
populations, the "gian t" glitchers and pulsa rs that only ex- 
hibit smaller glitches (JEspinoza et al.ll201ll ). We shall thus 
focus on the population o f "Vela-like" pul s ars th at show gi- 
ant glitches, as defined bv JEspinoza et al.l ((2011'), which are 
shown on the P — P diagram in figure [131 In particular 7 of 
these objects have multiple giant glitches, and one can thus 
derive an approximate waiting time between glitches. In fig- 
ure[T4]we plot the waiting time as a function of the spindown 
rate v. The data appears consistent with our theoretical ex- 
pectation that a pulsar that is spinning down faster will 
build up the critical lag on a shorter timescale and glitch 
more often. In particular we can see that most systems lie 
close to the Vela in the P—P diagram and also glitch roughly 
every few years, but the X-ray pulsar J0527-6910, which is 
spinning down approximately an order of magnitude faster 
glitches every few months, while the lower limits on slower 
pulsars indicate that they may glitch every decade. 

Naturally the situation will be complicated by the fact 
that these pulsars also exhibit smaller glitches, which may 
transfer part of the angular momentum before a giant glitch, 
and by t he fact that the critical lag will also depend, albeit 
weakly (|Pizzochero. Seveso fc Haskell l2011f ). on the mass 
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Figure 13. The location of th e giant glitching pulsars as iden- 
tified by lEspinoza et al.l l|201lh in the P — P diagram. This is 
the population of pulsars that show large steps in the frequency 
(Au Si 10~* Hz) and always exhibit an increase in the spindown 
rate after the glitch. 



and radius of the star. Given these limitation our model 
would, however, appear to be consistent with the observed 
inter-(giant)glitch waiting time. 

Unfortunately only the Vela is currently observed at 
short enough intervals to allow a fit for the short timescale 
transient terms in the relaxation, so more accurate tests are 
not currently possible for other pulsars. It would be of great 
interest if such observations for other giant glitchers were to 
become possible with the new generation of radio-telescopes 
such as LOFAR and the SKA. 



8 CONCLUSIONS 

We have presented a hydrodynamical two fluid model of 
pulsar glitches that can consistently model all phases of the 
glitch itself. Our model can be successfully applied to the 
giant glitches of the Vela pulsar, for which we can reproduce 
the approximate waiting time between glitches, the size of 
the glitch and the short term post-glitch relaxation. The 
main assumption is that a giant glitch will occur once the 
system exceeds the maximum lag that the pinning force in 
the crust can sustain. This naturally gives rise to a wait- 
ing time between glitches that depends on the pulsar spin- 
down rate (i.e. it is the time it takes the crust to spin down 
by the required amount) and to a maximum size for the 
glitch. Both these quantities depend only weakly on the mass 
and radius of the star ICrill fc Pizzochoro 2011; GriU 201l|; 
jPizzochero. Seveso fc Haskellll201ll ) and our model can, in 
fact, reproduce these features successfully also for the gen- 
eral population of "giant glitchers", i.e. the pulsars for which 
giant glitches such as those of the Vela have been observed 
(Espinoza et al. 2011j). 

In our model the coupling between the charged com- 
ponent (which we assume to be rigidly rotating) and the 
superfluid neutrons is given by the vortex-mediated mutual 
friction. Our results suggest that the mutual friction will be 
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Figure 14. We plot the approximate waiting time between 
glitches for the pulsars that have shown multiple giant glitches, as 
a function of the spin down rate. We also include two pulsars that 
have shown only one glitch but also have a long baseline for the 
observations and can thus provide us with an interesting lower 
limit on the waiting time. The data appears consistent with the 
notion that giant glitches can occur once a critical lag of approxi- 
mately Af2 = 10"'^ is reached. In fact the Vela like pulsars glitch 
every few years, but the X-ray pulsar J0527-6910, which is spin- 
ning down approximately an order of magnitude faster glitches 
every few months, while the lower limits on slower pulsars indi- 
cate that they may glitch every decade. In fact the data appears 
to be well described by a fit of the form y=A/x, as shown in 
the figure, with y the waiting time in seconds, x the frequency 
derivative and A = 1.082212 X IQ-^ Hz. 



model should relax the rigid rotation assumption for the 
charged component and include the effect of Ekman pump- 
ing. Further developments should also include more realistic 
models for the drag parameters in the star, as the density 
dependence of the coupling strength clearly has an impact 
on the amount of angular momentum that can be exchanged 
on different timescales. Truly quantitative results could then 
be obtained with the use of realistic equations of state to- 
gether with consistent estimates of t he pi nning force , such 
as those of (| Grill fc PizzocherdlJOlTI ) and (|Grillll201ll ). 

Note that we have assumed that a giant glitch only 
occurs when the maximum critical lag is reached. If unpin- 
ning could be triggered earlier, this could generate smaller 
glitches. In fact cellular automaton models have shown that 
the waiting time and size distributions of pulsar glitches can 
be successfully explained by vortex avalanche dynamics, re- 



Warszawski & Mclat 
Warszawski fc Melato; 



lated to random unpinning events 
2010l : iMelatos fc Warszawskil |2009| : 
2011f ). It would thus be of great interest to use our long-term 
hydrodynamical models, with realistic pinning forces, as a 
background for such cellular automaton models that model 
the short-term vortex dynamics. Such a model could then 
also be extended to model not only large pulsar glitches, 
but more generally pulsar timing noise, an issue that is of 
great importance for t he current efforts t o detects GWs with 
pulsar timing arrays (jHobbs et al.ll2010l ). 

Finally, the next generation of radio telescopes, such as 
LOFAR and the SKA, is likely to provide much more precise 
timing data for radio pulsars and is likely to set much more 
stringent constraints on the glitch rise time and short term 
relaxation, thus allowing us to test our models and probe 
the coupling between the interior superfluid and the crust 
of the NS with unprecedented precision. 



weak in the crust, possibly due to the fact that not all vor- 
tices are free, but rather that the strong pinning force gives 
rise to a situation in which most vortices are pinned and 
only a small fraction can 'creep' outwards. Only once the 
maximum unpinning lag is exceeded can the vortices move 
out freely; a process which can excite Kelvin oscillations and 
give rise to a strong drag and recoupling of the two compo- 
nents on a very short timescale, i.e. a glitch. The short term 
post-glitch relaxation of the Vela, on the other hand, sug- 
gests that the magnitude of the drag in the core of the NS is 
consistent with theoretical expectations for electron scatter- 
ing of magnetised vortex cores. Our model does not support 
the notion that, at least on short timescales, a significant 
number of vortices is pinned in the core (as could, for ex- 
ample, be the case if one has a type II superconductor and 
vortices cannot cross fluxtubes, effectively decoupling the 
core and the crust). A detailed analysis of the case in which 
the core consists of a type II superconductor will be a focus 
of future work in order to obtain more quantitative results 
and constraints on NS interior physics. Some vortices that 
cross the core may however be weakly pinned to the crust, 
and vortex repinning and creep (also in the core) may play 
a role on the longer timescales associated with the recovery. 
Another effect which will have an impact on the post- 
glitch recovery is the Ekman flow at the crust-core inter- 
face. This effect has been shown to be important in fitting 
the post glitch recovery of t he Vela and Crab pulsars by 
Ivan Evsden fc Melatod l|201Cl ) and future adaptations of our 
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